Project - Featurization, Model Selection & Tuning - AIML - Aishik Sengupta
Data Description :¶
The actual concrete compressive strength (MPa) for a given mixture under a specific age (days) was determined from laboratory. Data is in raw form (not scaled). The data has 8 quantitative input variables, and 1 quantitative output variable, and 1030 instances (observations).
Domain :¶
Cement manufacturing
Context :¶
Concrete is the most important material in civil engineering. The concrete compressive strength is a highly nonlinear function of age and ingredients. These ingredients include cement, blast furnace slag, fly ash, water, superplasticizer, coarse aggregate, and fine aggregate.
Attribute Information :¶
- Cement : measured in kg in a m3 mixture
- Blast : measured in kg in a m3 mixture
- Fly ash : measured in kg in a m3 mixture
- Water : measured in kg in a m3 mixture
- Superplasticizer : measured in kg in a m3 mixture
- Coarse Aggregate : measured in kg in a m3 mixture
- Fine Aggregate : measured in kg in a m3 mixture
- Age : day (1~365)
- Concrete compressive strength measured in MPa
Learning Outcomes:¶
- Exploratory Data Analysis
- Building ML models for regression
- Hyper parameter tuning
Objective :¶
Modelling of strength of high performance concrete using Machine Learning
Steps and tasks:¶
- Deliverable -1 (Exploratory data quality report reflecting the following) (20 marks)
- Univariate analysis (5 marks) - data types and description of the independent attributes which should include (name, meaning, range of values observed, central values (mean and median), standard deviation and quartiles, analysis of the body of distributions / tails, missing values, outliers
- Multivariate analysis (5 marks) - Bi-variate analysis between the predictor variables and between the predictor variables and target column. Comment on your findings in terms of their relationship and degree of relation if any. Presence of leverage points. Visualize the analysis using boxplots and pair plots, histograms or density curves. Select the most appropriate attributes
- Pick one strategy to address the presence outliers and missing values and perform necessary imputation (10 marks)
- Deliverable -2 (Feature Engineering techniques) (15 marks)
- Identify opportunities (if any) to create a composite feature, drop a feature etc. (5 marks)
- Decide on complexity of the model, should it be simple linear model in terms of parameters or would a quadratic or higher degree help (5 marks)
- Explore for gaussians. If data is likely to be a mix of gaussians, explore individual clusters and present your findings in terms of the independent attributes and their suitability to predict strength (5 marks)
- Deliverable -3 (create the model ) ( 15 marks)
- Obtain feature importance for the individual features and present your findings
- Deliverable -4 (Tuning the model) (20 marks)
- Algorithms that you think will be suitable for this project (5 marks)
- Techniques employed to squeeze that extra performance out of the model without making it overfit or underfit (5 marks)
- Model performance range at 95% confidence level (10 marks)
In [266]:
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
# various scaling algorithms
from sklearn.preprocessing import LabelEncoder
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import RobustScaler
from scipy.stats import zscore
from sklearn.impute import SimpleImputer
from sklearn.model_selection import KFold
from sklearn.model_selection import train_test_split # to split the dataset for training and testing
from sklearn.model_selection import cross_val_score # for performing cross validation of models
from sklearn.svm import SVC # Support vector machine model
from sklearn.utils import shuffle
from sklearn import metrics # to get various evaluation metrics
from sklearn.metrics import roc_auc_score # receiver operating curve score
from sklearn.metrics import accuracy_score # accuracy of prediction score
from sklearn.metrics import recall_score # recall score
from sklearn.metrics import precision_score # precision score
from sklearn.metrics import f1_score # f1 score
from sklearn.decomposition import PCA # performing principal components analysis
# Import Linear Regression machine learning library
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
from sklearn.tree import DecisionTreeRegressor
from sklearn.cluster import KMeans
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor, BaggingRegressor
from time import time
from scipy.stats import randint as sp_randint
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import RandomizedSearchCV
from sklearn.utils import resample
from sklearn.metrics import roc_auc_score
In [2]:
# read the dataset to a dataframe
cm_df = pd.read_csv('concrete.csv')
In [3]:
# taking a look at the first 10 rows of the dataframe, already seeing some missing values
cm_df.head(10)
Out[3]:
In [4]:
cm_df.shape # dimensions of the dataframe
Out[4]:
This is a regression problem containing all 8 predictor variables and a target variable.
In [5]:
# To check for null values
cm_df.isna().sum()
Out[5]:
There are no null values in the dataset provided to us
In [6]:
cm_df.info() # basic info such as datatype, value types
- All are numerical datatypes.
- Eight variables are hainvg 64 bit float datatype
- One variable is having 64 bit integer datatype rightly so since it counting the number of days
In [7]:
# Getting the 5 point summary
cm_df.describe().T
Out[7]:
- Many zero values can be observered for the predictors - ash, slag and superplastic
- Except slag, ash and age, other columns have almost same median and mean values.
In [8]:
#Checing how many rows have atleast one 0 value
(~cm_df.all(1)).sum()
Out[8]:
In [8]:
# To check for correlation
colormap = plt.cm.magma
plt.figure(figsize=(20,20))
sns.heatmap(cm_df.corr(), linewidths=0.1, vmax=1.0,
square=True, cmap=colormap, linecolor='white', annot=True)
Out[8]:
- Strength has a good correlation with cement, superplastic and age
- Water has a good correlation with fineagg and superplastic
In [9]:
cm_df.corr()
Out[9]:
In [10]:
#cm_df.median()
In [11]:
#To impute missing values
#impute = SimpleImputer( missing_values = 0, strategy='median')
#impute = impute.fit(df[:,0:18])
#cm_df.iloc[:,0:9] = impute.fit_transform(cm_df.iloc[:,0:9])
Exploratory Data Analysis¶
In [12]:
# Doing a boxplot on the whole data to check
sns.boxplot( data=cm_df, orient= "h" )
Out[12]:
In [13]:
plt.figure(figsize = (15,66))
cols = cm_df.columns.values
i=0 # column counter
j=1 # plot counter
k=1 # plot counter of each variable
while i < (len(cols) - 1):
if k == 1:
plt.subplot(9,2,j)
sns.distplot(cm_df[cols[i]])
j+=1
k+=1
plt.title(f'Distribution of {cols[i]}')
plt.xlabel(f'{cols[i]}')
plt.ylabel('Frequency')
else:
plt.subplot(9,2,j)
sns.boxplot(cm_df[cols[i]], color='yellowgreen')
j+=1
q1, q3 = np.percentile(cm_df[cols[i]],[25,75])
IQR = q3 - q1
plt.title(f'Boxplot of {cols[i]} \n \u03bc = {round(cm_df[cols[i]].mean(), 3)}, SE = {round(cm_df[cols[i]].std(),4)}, Median = {round(cm_df[cols[i]].median(),3)}, IQR = {round(IQR, 3)} ')
plt.xlabel(f'{cols[i]}')
plt.ylabel('Frequency')
i+=1
k=1
- Outliers can be observed in slag, water superplastic, fineagg ang age.
- A mix of gaussians can be observed in all predictors indicating that data was sourced from different places and consolidated here.
In [15]:
sns.pairplot(cm_df, diag_kind='kde')
Out[15]:
- Not much association can be gathered by looking at the pairplots.
- There is a seemingly linear relationship b/w cement and strength.
- Strength itself follows a near normal distribution.
Checking for outliers having a standard deviation of more than 3¶
In [14]:
from scipy import stats
import numpy as np
z = np.abs(stats.zscore(cm_df))
print(z)
In [15]:
print(np.where(z>3))
Capping to 3 SD for outliers in slag , water. fineagg and superplastic¶
In [17]:
for col in ['slag','water','fineagg','superplastic']:
percentiles = cm_df[col].quantile([0.01,0.99]).values
cm_df[col] = np.clip(cm_df[col], percentiles[0], percentiles[1])
#percentiles = cm_df['superplastic'].quantile([0.01,0.99]).values
#cm_df['superplastic'] = np.clip(cm_df['superplastic'], percentiles[0], percentiles[1])
Did not cap for age as the data is well beyond 3 SD and after doing some domain research, ageing is an important factor when it comes to strength. Its best left untouched at the moment.¶
In [18]:
#percentiless = cm_df['fineagg'].quantile([0.023,0.977]).values
#percentiless
In [19]:
#cm_df['slag'] = np.clip(cm_df['slag'], percentiless[0], percentiless[1])
In [20]:
#cm_df = cm_df[(z < 3).all(axis=1)]
In [21]:
plt.figure(figsize = (15,66))
cols = cm_df.columns.values
i=0 # column counter
j=1 # plot counter
k=1 # plot counter of each variable
while i < (len(cols) - 1):
if k == 1:
plt.subplot(9,2,j)
sns.distplot(cm_df[cols[i]])
j+=1
k+=1
plt.title(f'Distribution of {cols[i]}')
plt.xlabel(f'{cols[i]}')
plt.ylabel('Frequency')
else:
plt.subplot(9,2,j)
sns.boxplot(cm_df[cols[i]], color='yellowgreen')
j+=1
q1, q3 = np.percentile(cm_df[cols[i]],[25,75])
IQR = q3 - q1
plt.title(f'Boxplot of {cols[i]} \n \u03bc = {round(cm_df[cols[i]].mean(), 3)}, SE = {round(cm_df[cols[i]].std(),4)}, Median = {round(cm_df[cols[i]].median(),3)}, IQR = {round(IQR, 3)} ')
plt.xlabel(f'{cols[i]}')
plt.ylabel('Frequency')
i+=1
k=1
Now the outliers are gone
In [23]:
cm_df.corr()
Out[23]:
Composite Features¶
Added below new features after doing some domain research
In [24]:
cm_df['wc_ratio'] = cm_df['water'] / cm_df['cement']
cm_df['wb_ratio'] = cm_df['water'] / ( cm_df['water'] + cm_df['cement'] )
Not dropping columns at the moment
In [25]:
colormap = plt.cm.magma
plt.figure(figsize=(20,20))
sns.heatmap(cm_df.corr(), linewidths=0.1, vmax=1.0,
square=True, cmap=colormap, linecolor='white', annot=True)
Out[25]:
Thus we can see there is a very strong inverse correlation betweem the strength and the new features
In [26]:
cm_df.corr()
Out[26]:
In [27]:
sns.pairplot(cm_df, diag_kind='kde')
Out[27]:
The pairplots justify the addition of new features
Decide on the complexity of model¶
Order 1 complexity¶
Linear Models¶
In [28]:
X = cm_df.drop(['strength'], axis = 1)
y = cm_df['strength']
In [29]:
X_train, X_test, y_train, y_test = train_test_split(X , y ,test_size = 0.30 , random_state = 1234)
In [30]:
scale = StandardScaler() # Standard scaling
scale.fit(X_train.iloc[:,:])# fitting on training data so the data integrity in test is maintained
Out[30]:
In [31]:
X_train.iloc[:,:] = scale.transform(X_train.iloc[:,:])
X_test.iloc[:,:] = scale.transform(X_test.iloc[:,:])
In [32]:
X_train
Out[32]:
Linear Regression¶
In [33]:
reg_model = LinearRegression()
reg_model.fit(X_train, y_train)
for idx, col in enumerate(X_train.columns):
print("The coefficient for {} is {}".format(col, reg_model.coef_[idx]))
In [34]:
reg_model.coef_
Out[34]:
In [35]:
intercept = reg_model.intercept_
print("The intercept for our model is {}".format(intercept))
Ridge Regression¶
In [36]:
ridge = Ridge(alpha=.3)
ridge.fit(X_train,y_train)
print ("Ridge model:", (ridge.coef_))
Lasso Regression¶
In [37]:
lasso = Lasso(alpha=0.1)
lasso.fit(X_train,y_train)
print ("Lasso model:", (lasso.coef_))
Decision Tree Regression¶
In [38]:
dt_model = DecisionTreeRegressor( max_depth=10)
In [39]:
dt_model.fit(X_train, y_train)
Out[39]:
In [40]:
dt_model.score(X_test, y_test)
Out[40]:
In [41]:
dt_model.score(X_train, y_train)
Out[41]:
Comparing the scores¶
In [42]:
#Linear Regression
print(reg_model.score(X_train, y_train))
print(reg_model.score(X_test, y_test))
In [43]:
# Ridge Regression
print(ridge.score(X_train, y_train))
print(ridge.score(X_test, y_test))
In [44]:
# Lasso Regression
print(lasso.score(X_train, y_train))
print(lasso.score(X_test, y_test))
In [45]:
rfTree = RandomForestRegressor(n_estimators=50)
rfTree.fit(X_train,y_train)
print("rfTree on train data ", rfTree.score(X_train,y_train))
print("rfTree on test data ", rfTree.score(X_test, y_test))
Time to use Quadratic model to resolve non-linear interaction¶
In [46]:
#X_train = X_train.drop(['coarseagg','wc_ratio'], axis = 1 , inplace= True)
#X_test = X_test.drop(['coarseagg','wc_ratio'], axis = 1 , inplace= True)
In [47]:
from sklearn.preprocessing import PolynomialFeatures
poly = PolynomialFeatures(degree = 2, interaction_only=True)
In [48]:
X_train = poly.fit_transform(X_train)
X_test = poly.fit_transform(X_test)
In [49]:
X_train.shape
Out[49]:
In [50]:
X_test.shape
Out[50]:
In [51]:
#X_train, X_test, y_train, y_test = train_test_split(X_poly , y ,test_size = 0.30 , random_state = 1)
In [52]:
X_train
Out[52]:
In [53]:
X_train
Out[53]:
Linear Regression¶
In [54]:
reg_model = LinearRegression()
reg_model.fit(X_train, y_train)
Out[54]:
Ridge Regression¶
In [55]:
ridge = Ridge(alpha=.3)
ridge.fit(X_train,y_train)
print ("Ridge model:", (ridge.coef_))
Lasso Regression¶
In [56]:
lasso = Lasso(alpha=0.1)
lasso.fit(X_train,y_train)
print ("Lasso model:", (lasso.coef_))
Comparing the scores¶
In [57]:
#Linear Regression
print(reg_model.score(X_train, y_train))
print(reg_model.score(X_test, y_test))
In [58]:
# Ridge Regression
print(ridge.score(X_train, y_train))
print(ridge.score(X_test, y_test))
In [59]:
# Lasso Regression
print(lasso.score(X_train, y_train))
print(lasso.score(X_test, y_test))
In [60]:
rfTree = RandomForestRegressor(n_estimators=50)
rfTree.fit(X_train,y_train)
print("rfTree on train data ", rfTree.score(X_train,y_train))
print("rfTree on test data ", rfTree.score(X_test, y_test))
Decision Tree Regression¶
In [61]:
dt_model = DecisionTreeRegressor( max_depth=5)
In [62]:
dt_model.fit(X_train, y_train)
Out[62]:
In [63]:
dt_model.score(X_test, y_test)
Out[63]:
In [64]:
dt_model.score(X_train, y_train)
Out[64]:
In [65]:
X_train.shape
Out[65]:
Although Linear , Lasso & Ridge gave better results in order 2 complexity models, but order 1 is enough to capture the data. After fine tuning the models , we will resolve the overfit issue¶
Exploring K-Means Clusters based on mix of Gaussians¶
In [66]:
cm_df_z = cm_df.apply(zscore) # standardised score
In [223]:
# cluster analysis using Kmeans
cluster_range = range( 2, 15 )
cluster_errors = []
for num_clusters in cluster_range:
clusters = KMeans( num_clusters, n_init = 10 )
clusters.fit(cm_df_z)
labels = clusters.labels_
centroids = clusters.cluster_centers_
cluster_errors.append( clusters.inertia_ )
clusters_df = pd.DataFrame( { "num_clusters":cluster_range, "cluster_errors": cluster_errors } )
In [224]:
clusters_df # checking errors for various clusters created
Out[224]:
In [225]:
# Elbow plot to find optimal cluster
import matplotlib.pyplot as plt
plt.figure(figsize=(12,6))
plt.plot( clusters_df.num_clusters, clusters_df.cluster_errors, marker = "o" )
Out[225]:
Taking clusters 5 & 4 seems reasonable
Taking 5 clusters¶
In [227]:
kmeans = KMeans(n_clusters= 5)
kmeans.fit(cm_df_z)
Out[227]:
In [228]:
labels = kmeans.labels_
counts = np.bincount(labels[labels>=0])
print(counts)
In [234]:
## creating a new dataframe only for labels and converting it into categorical variable
cluster_labels = pd.DataFrame(kmeans.labels_ , columns = list(['labels']))
cluster_labels['labels'] = cluster_labels['labels'].astype('category')
cm_df_labeled = cm_df.join(cluster_labels)
cm_df_labeled.boxplot(by = 'labels', layout=(6,2), figsize=(20, 60))
Out[234]:
In [235]:
prediction=kmeans.predict(cm_df_z)
cm_df_z["GROUP"] = prediction # Creating a new column "GROUP" which will hold the cluster id of each record
cm_df_z_copy = cm_df_z.copy(deep = True) # Creating a mirror copy for later re-use instead of building repeatedly
Showing centroid data across various predictors for 5 cluster size¶
In [237]:
centroids = kmeans.cluster_centers_
centroids
Out[237]:
In [238]:
centroid_df = pd.DataFrame(centroids, columns = list(cm_df) )
centroid_df
Out[238]:
In [245]:
cm_df_z.boxplot(by = 'GROUP', layout=(4,3), figsize=(15, 25))
Out[245]:
Doing visual analysis on clusters for each predictor¶
In [242]:
var = 'cement'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[242]:
In [246]:
var = 'water'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[246]:
In [247]:
var = 'ash'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[247]:
In [248]:
var = 'slag'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[248]:
In [249]:
var = 'coarseagg'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[249]:
In [250]:
var = 'fineagg'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[250]:
In [251]:
var = 'age'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[251]:
In [252]:
var = 'superplastic'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[252]:
In [253]:
var = 'wc_ratio'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[253]:
In [254]:
var = 'wb_ratio'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z,hue='GROUP')
plot.set(ylim = (-3,3))
Out[254]:
Repeating same steps by taking 4 clusters¶
In [255]:
cm_df_z2 = cm_df.apply(zscore)
kmeans = KMeans(n_clusters= 4, random_state= 1)
kmeans.fit(cm_df_z2)
Out[255]:
In [256]:
labels = kmeans.labels_
counts = np.bincount(labels[labels>=0])
print(counts)
In [257]:
prediction=kmeans.predict(cm_df_z2)
cm_df_z2["GROUP"] = prediction # Creating a new column "GROUP" which will hold the cluster id of each record
cm_df_z2_copy = cm_df_z2.copy(deep = True) # Creating a mirror copy for later re-use instead of building repeatedly
In [258]:
centroids = kmeans.cluster_centers_
centroids
Out[258]:
In [259]:
centroid_df = pd.DataFrame(centroids, columns = list(cm_df) )
centroid_df
Out[259]:
In [260]:
cm_df_z2.boxplot(by = 'GROUP', layout=(4,3), figsize=(15, 25))
Out[260]:
In [261]:
var = 'cement'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[261]:
In [262]:
var = 'water'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[262]:
In [263]:
var = 'ash'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[263]:
In [264]:
var = 'slag'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[264]:
In [265]:
var = 'coarseagg'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[265]:
In [266]:
var = 'fineagg'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[266]:
In [267]:
var = 'age'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[267]:
In [268]:
var = 'superplastic'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[268]:
In [269]:
var = 'wc_ratio'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[269]:
In [270]:
var = 'wb_ratio'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=cm_df_z2,hue='GROUP')
plot.set(ylim = (-3,3))
Out[270]:
We can say that age, wc_ratio and wb_ratio are having linear relationship with strength , other than that Kmeans doesn't help much¶
Lets do some PCA¶
In [67]:
X = cm_df.drop(['strength'], axis = 1)
y = cm_df['strength']
In [68]:
scale2 = StandardScaler() # Standard scaling
scale2.fit(X.loc[:,:])# fitting on training data so the data integrity in test is maintained
X_scaled = scale2.transform(X.loc[:,:])
In [69]:
#should give an 18*18 matrix
covMatrix = np.cov(X_scaled,rowvar=False)
print(covMatrix)
In [71]:
#Performing PCA on all of its 10 predictors
pca = PCA(n_components=10)
pca.fit(X_scaled)
Out[71]:
In [74]:
plt.bar(range(1,11), pca.explained_variance_ratio_, alpha = 0.5, align='center', label='individual explained variance')
plt.step( range(1,11), np.cumsum(pca.explained_variance_ratio_), where='mid', label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc = 'best')
plt.show()
In [75]:
df_comp = pd.DataFrame(pca.components_,columns=list(X.columns.values))
In [77]:
df_comp
Out[77]:
In [76]:
plt.figure(figsize=(12,6))
sns.heatmap(df_comp,cmap='plasma')
Out[76]:
Building the Model to check for feature importances¶
In [142]:
X = cm_df.drop(['strength'], axis = 1)
y = cm_df['strength']
In [143]:
X_train, X_test, y_train, y_test = train_test_split(X , y ,test_size = 0.20 , random_state = 1)
In [144]:
X_train, X_val, y_train, y_val = train_test_split(X_train , y_train ,test_size = 0.25 , random_state = 1)
In [145]:
print(X_train.shape)
print(X_val.shape)
print(X_test.shape)
In [146]:
scale = StandardScaler() # Standard scaling
scale.fit(X_train.iloc[:,:])# fitting on training data so the data integrity in test is maintained
Out[146]:
In [147]:
X_train.iloc[:,:] = scale.transform(X_train.iloc[:,:])
X_val.iloc[:,:] = scale.transform(X_val.iloc[:,:])
X_test.iloc[:,:] = scale.transform(X_test.iloc[:,:])
In [148]:
gbmTree = GradientBoostingRegressor(n_estimators=50)
gbmTree.fit(X_train,y_train)
print("gbmTree on training" , gbmTree.score(X_train, y_train))
print("gbmTree on test data ",gbmTree.score(X_val,y_val))
print("gbmTree on test data ",gbmTree.score(X_test,y_test))
In [149]:
# View a list of the features and their importance scores
importances = gbmTree.feature_importances_
indices = np.argsort(importances)[::-1][:10]
a = cm_df.columns[:]
features= a.drop(['strength'],1)
#plot it
plt.figure(1)
plt.title('Feature Importances')
plt.barh(range(len(indices)), importances[indices], color='b', align='center')
plt.yticks(range(len(indices)), features[indices])
plt.xlabel('Relative Importance')
Out[149]:
In [150]:
rfTree = RandomForestRegressor(n_estimators=50)
rfTree.fit(X_train,y_train)
print("rfTree on train data ", rfTree.score(X_train,y_train))
print("rfTree on test data ", rfTree.score(X_val,y_val))
print("gbmTree on test data ",gbmTree.score(X_test,y_test))
In [151]:
# View a list of the features and their importance scores
importances = rfTree.feature_importances_
indices = np.argsort(importances)[::-1][:10]
a = cm_df.columns[:]
features= a.drop(['strength'],1)
#plot it
plt.figure(1)
plt.title('Feature Importances')
plt.barh(range(len(indices)), importances[indices], color='b', align='center')
plt.yticks(range(len(indices)), features[indices])
plt.xlabel('Relative Importance')
Out[151]:
Decided to drop below predictors based on multi-collinearity, Cluster analysis, PCA and feature imporatances¶
To use Train - Validation - Test Set of 3:1:1 ratio¶
In [152]:
X = cm_df.drop(['strength','coarseagg','ash','cement','wc_ratio'], axis = 1)
y = cm_df['strength']
In [153]:
X_train, X_test, y_train, y_test = train_test_split(X , y ,test_size = 0.20 , random_state = 1)
In [154]:
X_train, X_val, y_train, y_val = train_test_split(X_train , y_train ,test_size = 0.25 , random_state = 1)
In [155]:
print(X_train.shape)
print(X_val.shape)
print(X_test.shape)
In [156]:
scale = StandardScaler() # Standard scaling
scale.fit(X_train.iloc[:,:])# fitting on training data so the data integrity in test is maintained
Out[156]:
In [157]:
X_train.iloc[:,:] = scale.transform(X_train.iloc[:,:])
X_val.iloc[:,:] = scale.transform(X_val.iloc[:,:])
X_test.iloc[:,:] = scale.transform(X_test.iloc[:,:])
To use order 1 complexity regression models of SVR, Random Forrest, Bagging & Gradient Boost¶
In [158]:
svreg = SVR()
rfr = RandomForestRegressor()
gbr = GradientBoostingRegressor()
br = BaggingRegressor()
In [159]:
# Checking the hyperparameters present in the models for clarity
for clf, label in zip([svreg , rfr, gbr, br], ['svreg','rfr','gbr','br']):
print("model name: " , label)
print("\n model_hyperparameters \n" , clf.get_params() )
Time to use Grid Search & Random Search to fine tune models¶
SVR using GridSearchCV¶
In [36]:
svreg_gs_param_grid = {
'C' : [0.01, 0.1 , 1, 10,20, 30 , 50 , 100,200,400,500,1000],
'gamma' : ['auto','scale'],
'kernel' : ['poly','rbf']
}
In [37]:
# run grid search
grid_search = GridSearchCV(estimator=svreg, param_grid=svreg_gs_param_grid, cv=10)
#start = time()
grid_search.fit(X_train, y_train)
Out[37]:
In [38]:
#Getting best parameters
grid_search.best_params_
Out[38]:
In [39]:
print(grid_search.cv_results_['mean_test_score'])
In [40]:
grid_search.best_estimator_
Out[40]:
In [162]:
svregressor = SVR( C=200, gamma='auto')
svregressor.fit(X_train, y_train)
print("SVR on train data ", svregressor.score(X_train,y_train))
print("SVR on validation data ", svregressor.score(X_val,y_val))
print("SVR on test data ", svregressor.score(X_test,y_test))
Random Forrest Regression using GridSearchCV¶
In [44]:
rfr_gs_param_grid = {
'n_estimators' : [10, 50, 100, 200],
'max_depth': range(5,10),
'criterion': ['mse','mae'],
'min_samples_leaf' : range(1,4),
'max_features':['auto','sqrt']
}
In [45]:
# run grid search
grid_search = GridSearchCV(estimator=rfr, param_grid=rfr_gs_param_grid, cv=10)
#start = time()
grid_search.fit(X_train, y_train)
Out[45]:
In [46]:
# Getting the nest parameters
grid_search.best_params_
Out[46]:
In [163]:
rfTree = RandomForestRegressor(n_estimators=50, max_depth=11, max_features='auto',min_samples_leaf=1, criterion='mse')
rfTree.fit(X_train,y_train)
print("Random Forrest on train data ", rfTree.score(X_train,y_train))
print("Random Forrest on validation data ", rfTree.score(X_val,y_val))
print("Random Forrest on validation data ", rfTree.score(X_test,y_test))
Graddient Boost Regression Regression using GridSearchCV¶
In [63]:
gbr_gs_param_grid = {
'n_estimators' : [50, 100],
'max_depth': range(5,10),
'criterion': ['mse','mae'],
'min_samples_leaf' : range(1,4),
'max_features':['auto','sqrt'],
'learning_rate' : [0.001, 0.01, 0.05,0.1]
}
In [64]:
# run grid search
grid_search = GridSearchCV(estimator=gbr, param_grid=gbr_gs_param_grid, cv=10)
#start = time()
grid_search.fit(X_train, y_train)
Out[64]:
In [65]:
# Getting the best parameters
grid_search.best_params_
Out[65]:
In [164]:
gbmTree = GradientBoostingRegressor(n_estimators=100, criterion='mse', learning_rate=0.1, max_depth=6, max_features='sqrt',
min_samples_leaf=3)
gbmTree.fit(X_train,y_train)
print("Gradient Boost on training data" , gbmTree.score(X_train, y_train))
print("Gradient Boost on validation data ",gbmTree.score(X_val,y_val))
print("Gradient Boost on test data ",gbmTree.score(X_test,y_test))
Bagging Regression using GridSearchCV¶
In [26]:
br_gs_param_grid = {
'n_estimators' : [50, 100, 200, 400, 500, 1000],
'max_features':range(1,6),
'bootstrap':[True, False]
}
In [27]:
# run grid search
grid_search = GridSearchCV(estimator=br, param_grid=br_gs_param_grid, cv=10)
#start = time()
grid_search.fit(X_train, y_train)
In [70]:
# Getting the best parameters
grid_search.best_params_
Out[70]:
In [ ]:
print(grid_search.cv_results_['mean_test_score'])
In [165]:
bgcl = BaggingRegressor(n_estimators=500, bootstrap= True, max_features=5)
bgcl = bgcl.fit(X_train,y_train)
print("bgcl on train data ", bgcl.score(X_train,y_train))
print("bgcl on validation data ", bgcl.score(X_val,y_val))
print("bgcl on test data ", bgcl.score(X_test,y_test))
#print("out of bag score" , bgcl.oob_score_)
Bagging Regression using Random SearchCV¶
In [166]:
svreg = SVR()
rfr = RandomForestRegressor()
gbr = GradientBoostingRegressor()
br = BaggingRegressor()
In [167]:
br_rs_param_grid = {
'n_estimators' : [50, 100, 200, 400, 500, 1000],
'max_features':range(1,6),
'bootstrap':[True, False]
}
In [168]:
# run randomized search
samples = 60 # number of random samples
randomCV = RandomizedSearchCV(br, param_distributions=br_rs_param_grid, n_iter=samples) #default cv = 3
In [169]:
randomCV.fit(X_train, y_train)
print(randomCV.best_params_)
In [186]:
bgcl = BaggingRegressor(n_estimators=500, bootstrap= True, max_features=5)
bgcl = bgcl.fit(X_train,y_train)
print("bgcl on train data ", bgcl.score(X_train,y_train))
print("bgcl on test data ", bgcl.score(X_val,y_val))
print("bgcl on test data ", bgcl.score(X_test,y_test))
#print("out of bag score" , bgcl.oob_score_)
Gradient Boost Regression Regression using Random SearchCV¶
In [170]:
gbr_rs_param_grid = {
'n_estimators' : [50, 100, 200, 400, 500, 1000],
'max_depth': range(4,7),
'criterion': ['mse','mae'],
'min_samples_leaf' : sp_randint(1, 8),
'max_features':sp_randint(1,6),
'loss' : ['ls', 'lad', 'huber', 'quantile'],
'learning_rate' : [0.001, 0.01, 0.05,0.1, 0.2,0.3]
}
In [171]:
# run randomized search
samples = 1500 # number of random samples
randomCV = RandomizedSearchCV(gbr, param_distributions=gbr_rs_param_grid, n_iter=samples) #default cv = 3
In [172]:
randomCV.fit(X_train, y_train)
print(randomCV.best_params_)
In [174]:
gbmTree = GradientBoostingRegressor( criterion='mse', learning_rate = 0.05, loss= 'huber', max_depth= 5,
max_features= 2, min_samples_leaf= 7, n_estimators= 1000)
gbmTree.fit(X_train,y_train)
print("gbmTree on training" , gbmTree.score(X_train, y_train))
print("gbmTree on validation data ",gbmTree.score(X_val,y_val))
print("gbmTree on test data ",gbmTree.score(X_test,y_test))
Random Forrest Regression using Random SearchCV¶
In [180]:
rfr_rs_param_grid = {
'n_estimators' : [50, 100, 200, 400, 500, 1000],
'max_depth': range(5,10),
'criterion': ['mse','mae'],
'min_samples_leaf' : sp_randint(1, 5),
'max_features':sp_randint(1, 6)
}
In [181]:
# run randomized search
samples = 1000 # number of random samples
randomCV = RandomizedSearchCV(rfr, param_distributions=rfr_rs_param_grid, n_iter=samples) #default cv = 3
In [182]:
randomCV.fit(X_train, y_train)
print(randomCV.best_params_)
In [183]:
rfTree = RandomForestRegressor(n_estimators=500, max_depth=9, max_features=5 ,min_samples_leaf=1, criterion='mse')
rfTree.fit(X_train,y_train)
print("rfTree on train data ", rfTree.score(X_train,y_train))
print("rfTree on validation data ", rfTree.score(X_val,y_val))
print("rfTree on test data ", rfTree.score(X_test,y_test))
SVR using Random SearchCV¶
In [179]:
svreg_rs_param_grid = {
'C' : [0.01, 0.1 , 1, 10,20, 30 , 50 , 100, 200, 300, 500],
'gamma' : [0.01, 0.1,0.05,0.5,1,10,20, 'auto', 'scale']
}
In [60]:
# run randomized search
samples = 99 # number of random samples
randomCV = RandomizedSearchCV(svreg, param_distributions=svreg_rs_param_grid, n_iter=samples) #default cv = 3
In [61]:
randomCV.fit(X_train, y_train)
print(randomCV.best_params_)
In [185]:
svregressor = SVR( C=500, gamma=0.1)
svregressor.fit(X_train, y_train)
print("SVR on train data ", svregressor.score(X_train,y_train))
print("SVR on validation data ", svregressor.score(X_val,y_val))
print("SVR on test data ", svregressor.score(X_test,y_test))
In [252]:
y_test
Out[252]:
Gradient Boosting is the model to rely on as it gives better R2 results¶
Take a look at the feature importances¶
In [264]:
# View a list of the features and their importance scores
importances = gbmTree.feature_importances_
indices = np.argsort(importances)[::-1][:6]
a = cm_df.columns[:]
features= a.drop(['strength','coarseagg','ash','cement','wc_ratio'],1)
#plot it
plt.figure(1)
plt.title('Feature Importances')
plt.barh(range(len(indices)), importances[indices], color='b', align='center')
plt.yticks(range(len(indices)), features[indices])
plt.xlabel('Relative Importance')
Out[264]:
Bootstrap Sampling¶
In [253]:
Xdf = cm_df.drop(['coarseagg','ash','cement','wc_ratio'], axis = 1)
In [254]:
strength = Xdf['strength']
Xdf.drop(labels=['strength'], axis=1, inplace=True)
Xdf.insert(6,'strength',strength )
In [255]:
scale = StandardScaler() # Standard scaling
scale.fit(Xdf.iloc[:,:-1])# fitting on training data so the data integrity in test is maintained
Out[255]:
In [256]:
Xdf.iloc[:,:-1] = scale.transform(Xdf.iloc[:,:-1])
In [257]:
Xdf.head()
Out[257]:
In [258]:
val = Xdf.values
In [259]:
len(Xdf)
Out[259]:
In [260]:
val
Out[260]:
In [261]:
# configure bootstrap
n_iterations = 1000 # Number of bootstrap samples to create
n_size = int(len(Xdf) * 0.70) # picking only 50 % of the given data in every bootstrap sample
# run bootstrap
stats = list()
for i in range(n_iterations):
# prepare train and test sets
train = resample(val, n_samples=n_size) # Sampling with replacement
test = np.array([x for x in val if x.tolist() not in train.tolist()]) # picking rest of the data not considered in sample
# fit model
model = GradientBoostingRegressor( criterion='mse', learning_rate = 0.05, loss= 'huber', max_depth= 5,
max_features= 2, min_samples_leaf= 7, n_estimators= 1000)
model.fit(train[:,:-1], train[:,-1])
# evaluate model
#predictions = model.predict(test[:,:-1])
score = model.score(test[:,:-1], test[:,-1]) # caution, overall accuracy score can mislead when classes are imbalanced
print(score)
stats.append(score)
In [263]:
# plot scores
plt.hist(stats)
plt.show()
# confidence intervals
alpha = 0.95 # for 95% confidence
p = ((1.0-alpha)/2.0) * 100 # tail regions on right and left .25 on each side indicated by P value (border)
lower = max(0.0, np.percentile(stats, p))
p = (alpha+((1.0-alpha)/2.0)) * 100
upper = min(1.0, np.percentile(stats, p))
print('%.3f confidence interval %.3f%% and %.3f%%' % (alpha*100, lower*100, upper*100))
K-Fold Cross Validation¶
In [187]:
X = cm_df.drop(['strength','coarseagg','ash','cement','wc_ratio'], axis = 1)
y = cm_df['strength']
In [188]:
X.apply(zscore)
Out[188]:
Using whole dataset for K-Fold Validation¶
In [189]:
kfold = KFold(n_splits=10)
model = GradientBoostingRegressor( criterion='mse', learning_rate = 0.05, loss= 'huber', max_depth= 5,
max_features= 2, min_samples_leaf= 7, n_estimators= 1000)
results = cross_val_score(model, X, y, cv=kfold)
print(results)
print("Accuracy: %.3f%% (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
In [190]:
# To find interval at 95% confidence
std95 = results.std() * 1.96
min_interval = results.mean() - std95
max_interval = results.mean() + std95
In [195]:
print('We can inform the stakeholders that the model can be said to perform between {min:.3f}% & {max:.3f}% at 95% confidence'.
format(min = min_interval*100.0, max = max_interval*100.0))